Improved extremal optimization for the Ising spin glass 
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A version of the extremal optimization (EO) algorithm introduced by Boettcher and Percus is 
tested on 2D and 3D spin glasses with Gaussian disorder. EO preferentially flips spins that are 
locally "unfit"; the variant introduced here reduces the probability to flip previously selected spins. 
Relative to EO, this adaptive algorithm finds exact ground states with a speed-up of order 10 4 (10 2 ) 
for 16 2 - (8 3 -) spin samples. This speed-up increases rapidly with system size, making this heuristic 
a useful tool in the study of materials with quenched disorder. 
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Exploring the low temperature behavior of disordered 
materials, such as spin glasses and other random mag- 
nets Q , is quite challenging due to the very phenomena, 
glassy dynamics and multiple metastable states, that are 
important in such materials. Scaling arguments 
indicate that many properties of the glassy state, includ- 
ing the scaling of the energy of excitations and correla- 
tion functions, can be found by studying the ground state 
and its response to perturbations. Significant effort has 
been invested in identifying models whose ground states 
can be computed in time polynomial in the system size 
Q. Where no polynomial-time algorithm is known, ex- 
act and heuristic methods which take time exponential 
in system size are used. This enterprise is intimately 
connected with concepts developed in computer science, 
especially the distinction between P and NP-hard opti- 
mization problems 0|. 

The Ising spin glass (ISG) is a prototypical example of 
a disordered magnet. NP-hard problems such as the 3D 
ISG are, of course, particularly challenging. Exact meth- 
ods for the 3DISG with Gaussian bond weights can solve 
12 3 -spin samples with open boundary conditions Q. 
Such sizes have not proven to be sufficiently large to de- 
cide between alternate pictures for the low-temperature 
behavior. Heuristic genetic methods mix configurations 
and can therefore generate large scale "moves": such 
methods are used for samples with 14 3 spins for ± J cou- 
plings Jj. Heuristics with local moves generally have 
difficulty finding the exact ground state, due to the large 
barriers separating metastable states. Techniques such as 
flat histogram methods @ can partially lower free energy 
barriers between metastable states. 

In this Communication, I study a modified version of 
extremal optimization (EO) [10] . EO is a local search al- 
gorithm that preferentially flips spins with low "fitness". 
The version presented here, "jaded" extremal optimiza- 
tion (JEO) increases the fitness of a spin by an amount 
proportional to the number of times it has been flipped. 
The goal of this adjustment is to reduce the repetition 
in exploring paths in configuration space, so that more 
possibilities can be quickly explored. Empirically, this 
simple change dramatically increases the effectiveness of 
the EO algorithm for finding ground states of two- and 
three-dimensional spin glass samples. As exact ground 
states are needed for studies of excitations and scaling, 
the algorithm is, for the most part, stringently tested by 



demanding that it find the ground states computed by 
exact methods. Both EO and JEO take time exponen- 
tial in the system size to find the exact ground state, 
but the rate of growth is slower for JEO. Though JEO 
introduces an extra parameter, large improvements are 
achieved with only modest tuning. 



I. EXTREMAL OPTIMIZATION AND 
EXTENDED ALGORITHM 

A principle motivation for applying EO is to explore 
the energy landscape near the trial configuration by un- 
conditionally modifying "unfit" variables. Preferentially 
(but not exclusively) changing variables with low fitness 
tends to raise the expected fitness while maintaining large 
fluctuations. The algorithm differs some from traditional 
Monte Carlo algorithms that conditionally select vari- 
ables according to the expected improvement. In EO, the 
potential moves are selected according to their rank by 
fitness, rather than a Boltzmann distribution by weight. 

A correspondence can be defined between fitness and 
the Hamiltonian for the Ising spin glass The 
Hamiltonian for spins Si, indexed by position i, in a d- 
dimensional ISG of linear size L is 



H 



(ij) 



(1) 



where Jij are random bond strengths each chosen with 



probability P{Ji : 



e~ J ij/ 2 /y // 2n for nearest neighbor 
spins with 1 < i, j < N — L d . When d = 2, algorithms 
with running times polynomial in N are available [TlJ to 
find the ground state. When d > 3, finding the ground 
state energy is NP-hard, so that finding ground states 
for the worst-case choice of Jij is expected to take time 
exponential in N. In the context of EO, one choice for 
the fitness variable Xi for a spin variable Sj is 
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where Ui are the set of unsatisfied bonds (siJijSj < 0) 
containing Si. (Allowing for site-dependent constant 
shifts A- 1 — > A" + m as in Ref. 0] did not affect the 
comparisons here.) The configuration energy is related 
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to the fitness by H — — | J^i A° + I -A? I • Any increase 
in the fitness decreases the total energy. 

Given the fitness variables A^, there are a variety of 
strategies one could employ to attempt to improve the 
total fitness. The simplest version of EO takes "greedy" 
steps: the algorithm repeatedly flips the least fit vari- 
able until a static state is achieved. The greedy method 
converges quite rapidly, but in a spin glass the conver- 
gence is to a local minimum that is generally quite far 
from the optimal solution, both in configuration of the 
{si} and often in energy per degree of freedom H/N. 
Similar greedy approaches for decision problems such as 
SAT, which seeks truth assignments for Boolean formula 
so that all clauses contain a true value, can be quite suc- 
cessful for given ensembles of problems 

An improved method, r-EO sorts the spins by 

Xi and chooses the mth spin in the list with probabil- 
ity proportional to m" r . This favors the choice of spins 
with low fitness, but allows for the occasional choice of 
sites with very high fitness. Fluctuations arising from 
the stochastic choice among spins with low fitness and 
the ranking of spins by the total weight of broken bonds, 
rather than energy improvement, allow the search to es- 
cape metastable states. It is argued that for large 
systems, the optimal choice of t approaches t = 1. 

The extension considered in this paper (JEO) adjusts 
the fitness by an amount proportional to the number of 
times ki that a site i has been previously chosen, that is, 

Xi = Af = A° + Tki, (3) 

where T is a site-independent "aging" parameter. The 
variables are sorted by Af and then selected by rank as 
in t-EO. The t-EO algorithm corresponds to the choice 
T = 0. Setting T ^ reduces the probability of selecting 
moves that have been flipped many times before. For con- 
figurations near (or in) the ground state, it is favorable 
for some spins to have low fitness, in order that a number 
of other spins can maximize their fitness. When T = 0, 
these spins, which are actually in their ground state ori- 
entation relative to the other spins, will be flipped in fu- 
tility. Shifting the Xi during the algorithm also breaks the 
finite set of offsets between fitnesses of distinct spins that 
exist at r = (due to the finite number of bond configu- 
rations at each site). This adaptive scheme has similari- 
ties to a variety of methods for solving problems such as 
SAT (satisfiability of sets of logical constraints) that dis- 
favor repeated selection of the same move, such as Nov- 
elty 0] and variants of WALKS AT and GSAT (HGU- 
In contrast with these other schemes, the selection pro- 
cess in JEO is combined with the power law distribution 
for selecting ranked moves. Spin glasses with continuous 
disorder differ from SAT problems as they have less local 
degeneracy but also possess a global up-down symmetry, 
so that distinct methods may be appropriate. 

In order to select spins quickly, I used the approximate 
selection method described in Ref. [Tlj . The spins are 
stored in a heap structure 0] according to their current 
fitness. This structure is a tree that is relatively cheap to 



maintain (O(logTV) total cost to select a spin and update 
the tree). Each spin has a parent (except for the root) 
and at most two children. Each child is more fit than 
its parent and the root of the tree contains the least fit 
spin. This structure does not guarantee any other inter- 
level sorting, so that a spin i that is deeper in the tree 
than, but not a direct descendant of, a given spin i' , may 
have a lower fitness. The heap structure does maintain a 
useful approximate sorting, though. To select a spin to 
flip, a level £ is selected with probability proportional to 
2~( T ~ 1 ^ and then a random spin within level i is chosen. 
The spin at this site is then inverted. The fitness of the 
neighboring spins is adjusted and the heap is updated 
using standard methods [17| . 

EO does not take advantage of the special structure 
of the 2D problem: it is not necessary or even expected 
that it will find the solution in time polynomial in the sys- 
tem size. Polynomial-time solvable problems have been 
used to study algorithms, for example, for hard mean- 
field problems yj§. For some classes of problems , he uris- 
tics can find solutions in polynomial time fl3l UJ. In 
the 2DISG, large low-energy excitations may make local 
algorithms especially inefficient. 



II. PERFORMANCE OF THE ALGORITHM 

In this section, I compare the performance of the ex- 
tended EO algorithm, JEO, against r-EO as applied Ising 
spin glasses with Gaussian disorder. When feasible, com- 
parisons with ground states found using exact methods 
provide a precise and direct test for convergence. 

Two-dimensional spin glass. The 2DISG models are on 
a square lattice with L 2 spins and open boundary con- 
ditions. To determine the 2D ground state, each sample 
is mapped 0] to a general weighted matching problem. 
The matching problem for a graph is to find a set of 
edges with minimal total weight such that each vertex 
belongs to exactly one edge. The weighted graph for a 
2DISG sample has edges dual to the lattice bonds, with 
weight \Jij\ for an edge that crosses a bond with weight 
Jij , and extra edges of weight zero that ensure that the 
frustration of each plaquette is maintained: unfrustrated 
(frustrated) plaquettes give an even (odd) number of the 
bonds dual to the edges of the plaquette in the match- 
ing. To find the minimum weight matching and hence 
the ground state energy for a 2DISG sample, I used the 
Blossom IV algorithm developed by Cook and Rohe [2(i| . 

The exact ground state energy of each 2DISG sample 
was input to the r-EO and JEO codes. When the heuris- 
tic codes found this energy, the codes terminated. The 
primary results from these computations were the distri- 
butions of the running times, measured in number of spin 
flips, to find the true ground state. The time to solution 
is a function of both the seed used to generate the sample 
and an independent "algorithm seed" used to generate the 
random initial configuration and to select spin flips. In a 
given sample, the distribution of times to find a ground 
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Figure 1: Plot of t m , the sample mean of the median time 
to find the ground state , measured in spin flips, using r-EO 
(squares) and JEO (circles), for the 2DISG with optimal r 
and, for JEO, T. The triangles indicate the same measure of 
time to find the ground state energy to within 1% accuracy. 
The line shows, for comparison, a running time exponential 
in L, t m — 15 ■ 2 L , consistent with the results for JEO. The 
uncertainties are comparable to the symbol size. 



state was roughly Poissonian. This suggests that restart- 
ing the algorithm with different initial configurations or 
seeds for selecting flips does not significantly decrease the 
mean running time. This conclusion was consistent with 
empirical trials of restarting the algorithm: the algorithm 
does not get stuck in history dependent traps. Given a 
sample k, the median t„ of the running time was esti- 
mated from the solution time for 100 algorithm seeds. 
The results reported here are for t m , the sample mean 
of t™,. The r = data is in agreement with previously 
results for r-EO, with t m minimal at r w 1.5. 

The results for the mean solution time t m for optimal 
r and T are summarized in Fig. As suggested by the 
data plotted in Fig.0 t m is not very sensitive to the exact 
choice of parameters, as long as r is in the range 1.5 < 
r < 2.5 and the optimal Y (on the order of 10~ 3 to 10 _1 ) 
is found to within a factor of about 2, for the sizes studied 
here. The best running times for t-EO grow much more 
rapidly than those for JEO. For L = 16, JEO is of the 
order 10 4 times faster than t-EO. Extrapolation suggests 
that the advantage of JEO increases significantly with L. 
For comparison, an exponential dependence t m = 15 • 2 L 
is shown in Fig. This function does a good job of 
describing the JEO data for L = 4 through L = 32. 
In separate runs, for comparison, the heuristic algorithm 
was terminated when the energy was within 1% of the 
exact ground sate energy. These approximate solutions 
were found much more rapidly than exact solutions (~ 
10 5 times faster for L = 32). 

Three-dimensional spin glass. A similar comparison 
was carried out for 3DISG samples with Gaussian disor- 
der. The L 3 spins in the 3DISG samples lie on a cubic 
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Figure 2: Plot of t m for 2DISG samples of size L = 8, for F 
ranging from T = (i.e., r-EO) through T = 0.5, as a function 
of the power law for rank selection, r. For clarity, the error 
bars, which are of order 10% of the values for all points, are 
not shown. The solid lines are added only to group the points. 
Choosing F ~ 0.1 and r w 2.0 minimizes the run time. 



lattice with periodic boundary conditions. For 3DISG 
samples of size up to 6 3 , the spin glass server at the Uni- 
versity of Koln [2]J (which applies branch- and-cut @) 
was used to generate exact solutions. The termination 
condition of the algorithm was modified, as exact ground 
states for the larger samples were not readily available. 
All samples were simulated in parallel with n = 10 algo- 
rithm seeds. When the minimal record energy for eight 
(8) of the samples were identical, the algorithm was ter- 
minated. This criterion produced configurations equal 
to the exact solutions for all L — 4,6 samples (45 at 
each size). This suggests that true ground states were 
found with a high probability for L = 8 and possibly 
also L = 10. The summary results are plotted in Fig. 
Given the termination criterion, JEO was of the order of 
10 2 times faster than r-EO in converging to a potential 
solution for L = 8 samples. Very roughly, L = 6 samples 
were solved in ss 10 s on average both on the Koln spin 
glass server (a 400 MHz Sun Ultra) and using JEO (on 
a 1 GHz Intel P5). Further studies would be needed to 
provide better estimates of the confidence in the ground 
states and how to improve such confidence. 



III. DISCUSSION 

JEO extends the extremal optimization algorithm of 
Boettcher and Percus by adaptively reducing the fre- 
quency of flipping previously selected spins. As a local 
move can lead to avalanche-like behavior, due to induced 
changes in the fitness of neighbors, this modification also 
reduces the frequency of flipping larger domains. This 
extension of EO does add a parameter, the aging pa- 
rameter r. However, a near-optimal value for T for each 
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Figure 3: Plot of the sample average of the median running 
times for r-EO (squares) and JEO (circles) for the Gaussian 
Ising spin glass on a cubic lattice. The algorithm terminated 
when 8 of the minimal record energies agreed among 10 par- 
allel samples. The parameter r was fixed for JEO at a near- 
optimal r = 1.7 and near-optimal values of T = 0.1,0.1,0.05 
for L — 4,6, 8, respectively, were used. The gain for JEO over 
r-EO is approximately a factor of 100 at L = 8. The line 
shows t m = 0.05 • 2 3 ' 4 L , for a rough comparison. 

problem type at a given size can be found quickly and 
less tuning of the parameter r is required than for r-EO. 

One possible avenue of exploration is to check whether 
avalanche regions correspond to important domains or 



excitations in the sample. Possible modifications of JEO 
include using a selection distribution with sharp cutoffs 
[2^ . rather than power-law distributions. Other schemes 
for reducing the fitness of frequently repeated moves 
could be considered, such as modifying the fitness using 
non-linear functions of the number of flips at a spin. 

Regardless of the exact details of the role of domains 
and possible improvements, empirical testing shows that 
the aging of the spins during state-space exploration 
greatly reduces the time for EO to find the ground state 
of the ISG in two and three dimensions. Though the 
2D model was used to make a precise comparison with 
exact results, the exponential equilibration times for the 
2DISG using extremal optimization are consistent with 
those that would be seen for an NP-hard optimization 
problem with a similar local solution strategy. It may be 
useful to use an algorithm like JEO to locally improve 
the configurations formed by whole sample crossover in 
genetic algorithms [23] . As exact solutions for small sam- 
ples can be found with confidence in a relatively small 
number of steps, in machine time very similar to that 
for branch- and-cut, this simple algorithm also provides a 
very convenient way to study small 3D samples. 
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